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We investigate the properties of the molecular quantum dot (Holstein-Anderson) model using 
nnmerical and analytical techniques. Path integral Monte Carlo simulations for the cumulants of the 
distribution function of the phonon coordinate reveal that at intermediate temperatures the effective 
potential for the oscillator exhibits two minima rather than a single one, which can be understood as 
a signature of a bistability effect. A straightforward adiabatic approximation turns out to adequately 
describe the properties of the system in this regime. Upon lowering the temperature the two potential 
energy minima of the oscillator merge to a single one at the equilibrium position of the uncoupled 
system. Using the parallels to the X-ray edge problem in metals we derive the oscillator partition 
function. It turns out to be identical to that of the Kondo model, which is known to possess a 
universal low energy fixed point characterized by a single parameter - the Kondo temperature Tk- 
We derive an analogon of Tk for the molecular quantum dot model, present numerical evidence 
pointing towards the appearance of the Kondo physics and discuss experimental implications of the 
discovered phenomena. 

PACS numbers: 73.21.La, 05.10.Ln, 71.38.-k, 67.85.-d, 67.85.Pq 


In view of the recent progress in the field of microelec¬ 
tronic fabrication, which produces ever smaller electronic 
circuitry elements, it is reasonable to assume that the ba¬ 
sic building blocks of the future nanoelectronics would be 
individual molecules.® Contrary to the solid-state based 
systems their internal degrees of freedom play a princi¬ 
pal role. The m ost i mportant ones are the vibrational de¬ 
grees of freedom.®^ Although it is possible to model their 
effects with the help of a rather simple model - the molec¬ 
ular quantum dot (sometimes also referred to as Holstein- 
Anderson model), its properties are still not understood 
in full detail [1] . One of the reasons is that the problem in 
general is not exactly solvable and many of the interest¬ 
ing regimes are not accessible analytically. In particular, 
some time ago it was predicted that when the electron- 
phonon coupling is sufficiently strong, such a molecular 
dot might possess a bistability regime.®® A subsequent 
numerical analysis of systems under nonequilibirum con¬ 
ditions (with a finite voltage bias applied across the dot) 
has revealed some signatures of this phenomenon. How¬ 
ever, so far the numerics were not able to supply conclu¬ 
sive evidence about the lifetime of the s yst em in different 
conformational states of the molecule.®^ On the other 
hand, there are several arguments against a bistable be¬ 
haviour of such systems at low energies The pur¬ 
pose of this paper is to reconsider the problem, trying to 
settle the open issues outlined above for systems in equi¬ 
librium. By doing this we have made a twofold progress. 
Firstly, we report path integral Monte Carlo (PIMC) sim¬ 
ulations for the coordinate distribution functions of the 
localized vibrational degree of freedom, which are espe¬ 


cially convenient for measurements in future experimen¬ 
tal realizations of the model with the help of ultracold 
gas systems. Furthermore, we show a natural extension 
of the adiabatic approximation, which is valid in the low 
energy sector. There, surprisingly, the system undergoes 
a crossover into a regime which closely resembles the low 
energy limit of the Kondo model and which is determined 
by a single energy scale. 

The model for the molecular quantum dot possesses a 
single electronic level with energy ej and one vibrational 
mode, see Fig. (we consider a spinless system and use 
units in which H = e = 1). The corresponding local 
phonon is just a harmonic oscillator with mass m and 
frequency ujq. It is linearly coupled with strength A to 
the molecular electronic state which, in turn, is coupled 
via a tunnelling amplitude 7 to the fermionic continuum 
in a metallic electrode kept at chemical potential /i = 0 
at inverse temperature (3. The Hamiltonian comprising 
the aforementioned features reads 

iJ = Flel + i?ph + i?int , (1) 

He\ = Cdd) d + '^\ekc\ck + {a c\d + lA.c^ , 
k 

~ 2 m {d^d — 1 / 2 ). 

Cfc and d are the fermionic annihilation operators corre¬ 
sponding to the electronic states of the metallic electrode 
and of the molecule, respectively. P and Q are the mo¬ 
mentum and position operator of the phonon. 

Our main goal is to derive a numerical procedure to ac¬ 
cess the probability distribution of the expectation value 
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FIG. 1. (Color online) Schematic representation of the system 
described by the Hamiltonian 0 . 
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of the phonon position operator for arbitrary parame¬ 
ters, e. g. various coupling strength A and temperatures 
T, to search for signatures of bistability. In order to 
achieve this goal we first produce an effective action for 
the phonon by integrating out the electrode degrees of 
freedom out of 0 using standard functional integration 
techniqne. The resulting action is given by 


S'= S'o[Q]-h / (irdr'd'^{ t)Dq^{t — T')d{T') 

Jo 

rP 

-A / dr Q{t) d^T)d{T), (2) 

Jo 

where So[Q] is the action for a free harmonic oscillator 
and where by abuse of notation d^(r) and d(r) denote 
adjoint variables. Dq{t) is the Matsubara Green’s func¬ 
tion (GF) of the dot electron level, which can be found 
by elementary means. We concentrate on the wide flat 
band model, in which the conductance band of the elec¬ 
trodes is infinitely wide and has a constant density of 
states po-^Then one obtains 

TIq (Ain) “ ^Aln Cd IF Sgn rin (3) 


in energy representation with = 7r(2n l)/j3 and 

F = TrpoIdP being the inverse lifetime of the electron 
on the dot. Further integrating out the dot’s electronic 
degree of freedom yields for the partition function: 


Z = J 25Q e“‘5o[Q]-5 indet[p5„„£) 


This expression cannot be evaluated analytically. It 
is, however, amenable to numerical treatment and is the 
basis for our PIMC simulations. We follow the procedure 
outlined in [HHni and adapt it for measurement of the 
irreducible cumulants of the type q = (Q) and = 
{{Q— (Q))") for n > 1 . We use the notation for 

the corresponding distributions, which are the principal 
quantities of interest. 

The notation f DQ within the imaginary-time path 
integral Q symbolically refers to the continuum limit 

fc=l 


FIG. 2. Dependence of the phonon displacement probability 
distribution on electron-phonon coupling strength A obtained 
via PIMG simulations for PV = wq/F = 1. A increases from 
zero (curve with the highest maximum) to A = 6 (curve with 
the outmost maxima) in unit steps. In all simulations muJo = 
1 so that go and A are measured in units of y/VJrmjjQ and 
■y/FmojQ, respectively. Dots represent the numerical data and 
curves are fits of a symmetric and normalized superposition 
of two Gaussians. In all simulations Cd = 0. 


where K is the number of steps the imaginary time span 
of interest was devided into and Ar is the corresponding 
step width. Qk indicates the phonon position at step k. 
Af denotes a normalization factor. It drops out of any 
observable since they all are given by ratios of limits of 
the type 0. For numerical purposes, the aforementioned 
limit is not performed, leaving the discrete version of 
0 to be a high yet finite dimensional integral which 
may readily be evaluated by the means of Monte Garlo 
methods. The numerically evaluated mean of the phonon 
displacement Q to some power n for example reads 


K 




E 


(Qfc-Qfc-l) , 2,^2 

—(aTP —+‘^oQk 


k=l 


( 6 ) 


X det ^ 0 ,n—m )] ; 


where is the Fourier transform of the thermal 

analogon Gj{ t) of the time-ordered GF. The latter is 
connected to the Matsubara GF Dq{t) by Keldysh rota¬ 
tion, resulting for our model in 


G^(r)=0(T)G>(r) + 0(-r)G<(r) 
duj 2 Fe-‘^^ 


Go>(r) = 


Go<(r) = 


Loo (w - Cd)^ + r2 
/■°° duj 2Fe-‘^^ 
'_oo - ed)^ + r2 


[/(^)-l] , (7) 

/(w) ■ 


Here, /(w) denotes the Fermi distribution function of the 
electrode. Note that Gq (r) is only finite for —/3 < t < 0 
and Gg (r) for 0 < r < /3. The Heaviside step function 
must be 0(0) = 1/2 at the origin. 
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Applying Monte Carlo method to § in order to es¬ 
timate it simply means stochastically sampling its inte¬ 
grand ~ i. e. randomly travelling phase space and eval¬ 
uating the integrand at the points visited. If samples 
were drawn uniformly (which for an infinite interval is 
rather a challenge), numerical effort would scale expo¬ 
nentially with the dimension of the space one wishes to 
sample. We employ an importance sampling procedure 
instead, which uses the fact that the physical distribution 
functions ( 9 *-"^) actually have a rather small support 
- interactions and Boltzmann factors reduce the rele¬ 
vant phase space volume significantly. Thus samples may 
not be drawn uniformly but according to a distribution 
-f’Mc({Qfe})) such that only the dominant areas of the 
phase space are sampled with high accuracy. The non- 
uniform exploration of phase space must be countered by 
accumulating the ratio of the integrand and the proba¬ 
bility density it is drawn from instead of the plain inte¬ 
grand. This strategy is applied optimally if TMc({Qfc}) 
is identical to the integrand of the integral to be esti¬ 
mated. In our case this optimum can be reached since 
due to the thermal type of our problem the integrand 
is non-negative and normalizable. However, this cannot 
be done in case of a system in non-equilibrium, when, 
for instance, the quantum dot is coupled to two metallic 
electrodes with different chemical potentials. In order to 
produce samples distributed according to the integ rand 
of (|^ we employ the Metropolis-Hastings algorithm. h^ l ^^ l 
That is, we construct a Markov chain of paths {Qk} 
whose equilibrium distribution is equal to the one we are 
looking for. Thus, once equilibrated, the Markov walk 
produces a manifold of trajectories, since any step of it 
represents an entire phonon path. 

Both the systematic error due to the finite number 
of discretization steps and the statistical error due to 
stochastic nature of simulations can be made arbitrar¬ 
ily small by increasing the discretization resolution and 
number of samples, respectively. In addition, Markov 
chain methods are a source of yet another error type. 
The obtained random samples are not statistically in¬ 
dependent since they are not drawn from a distribution 
but constructed from the Markov walk. This walk paces 
configuration space, rendering its consecutive steps cor¬ 
related. That implies a twofold difficulty: first of all, 
variances are underestimated if calculated in the same 
way as for independent numbers and, secondly, there 
is no means of assessing the convergence of individual 
walks. Only if two points are separated by more then 
the autocorrelation time, they can be considered being 
statistically independent. In order to control that we 
use the procedure proposed by Flyvbjerg,!^ according to 
which samples are paired, pair averages are treated as in¬ 
dependent samples and the statistical error is estimated. 
Subsequently, the just obtained averages are paired once 
again and the error is estimated anew, and so forth. The 
iteration is stopped, once the estimators form a plateau. 
In this way we obtain very reliable numerical data, the 
overall error of which is extremely small. 


Fig. shows the distribution of the average value (Q) 
for increasing electron-phonon coupling. For A = 0 we 
deal with a simple harmonic oscillator and Pi (q) is given 
by the probability distribution of the groundstate wave 
function. For growing interaction strength one observes 
a splitting of the central maximum into two resulting in a 
bimodal distribution. This can be interpreted as a forma¬ 
tion of a double dip in the oscillator’s effective potential, 
the latter in close analogy to the uncoupled case, where 
one would expect Pi{q) ~ exp[—/3I/(g)] in the semiclas- 
sical limit, with V{q) denoting the oscillator’s potential 
energy. This is closely related to the bistability effect dis¬ 
cussed earlier.l^ An approximation used in these works 
assumes the phonon degrees of freedom to be much slower 
than the electronic ones, e. g. when the typical timescale 
for the electron dynamics on the dot ^ F“^ is much 
smaller than the phonon oscillation period ^ F”^ 
one can use the adiabatic, or Born-Oppenheimer approx¬ 
imation. Then the field Q{t) in the A-term in 0 can be 
taken to be static: — Ilm) = SnmPQ with constant 

Q. The resulting functional integral can then be evalu¬ 
ated by virtue of the saddle point approximation, which 
is performed in , yielding in the case of strong coupling 
A^ 3> ttiujqT for Pi{q) two 5-shaped maxima at the po¬ 
sitions Q± = ±(7max where qmax = X/2mujQ. It turns 
out that a careful evaluation of the adiabatic approxi¬ 
mation beyond the saddle point is indeed able to recover 
the actual shape of the distribution function obtained 
numerically with high degree of accuracy. Moreover, the 
prediction for Q± is perfectly reproduced by our simula¬ 
tions. In general the central dip in Pi{q) at the original 
equilibrium position of the uncoupled oscillator q = 0 
might or might not touch the axis of ordinates. We call 
the former case perfect bistability while in the latter sit¬ 
uation we are dealing with the spurious bistability. 

While in the uncoupled case the distribution 
of higher moments for n > 1 can trivially be found from 
Pi{q) of q = {Q) via the relation 

it is not possible any more in the case A yf 0 because 
in a fully interacting system a multitude of irreducible 
multi-particle correlations emerge and the single-particle 
picture, which is essential for the derivation of Q breaks 
down. This is demonstrated in Fig. where we have 
plotted the distribution P 2 {q^^^) of the second cumulant 
for an interacting system as well as the prediction 

Upon further lowering the temperature one observes a 
rather fast collapse of the bimodal distribution function 
to one with a single maximum at Q = 0. This kind of 
crossover occurs at a very low critical temperature, which 
strongly depends on A and can be seen in PIMC results 
for two different values of the coupling, A = I and 2.5. 
Fig. a shows the data for intermediate coupling strength 
A = 1. Here the distribution at /3F = 2 is bimodal due 
to the bistability effect, while Pi{q) exhibits one unique 
maximum for all investigated temperatures above or be¬ 
low /3F = 2. The collapse of the spurious bistability at 
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FIG. 3. Distribution function of the variance for different 
electron-phonon strengths: A = 0 (the corresponding data 
points perfectly coincide with the prediction of ([^, plotted as 
a curve), and A = 1, 2, 3,4, 5, 6 (the maxima of the respective 
data sets lie at larger for growing A). 



FIG. 4. Phonon displacement probability distribution for 
A = 1 and uiq/T = 1 for different inverse temperatures /3. 



low temperature can be even more clearly observed in a 
strongly coupled system with A = 2.5, see Fig.[^ It turns 
out that under these conditions the bistability starts to 
emerge already at j3T « 0.2 and becomes almost per¬ 
fect at PT « 2 (in a sense that then P(0) « 0, which 
means that the probability to find the phonon localized 
in the original equilibrium position of the uncoupled os¬ 
cillator vanishes completely), persisting this way all the 
way down to temperatures f5T « 15. After that the bista- 
bilty starts to degrade as demonstrated by the curves for 
/3r « 20 and PV « 30 , which no longer touch the ordi¬ 
nate axis at (3 = 0. 

This interesting behaviour leaves open a number of im¬ 
portant questions. The first and most important of them 
regards the dependence of the lower crossover tempera¬ 
ture Tc, below which T’i(g) again shows one unique max¬ 
imum, on the electron-phonon coupling strength. The 
adiabatic approximation is not able to estimate T^. In 
what follows we improve it by taking care of the low- 
energy fluctuations of the Q{t) field in the functional 
integral Q. In order to make progress it is helpful to 
return to the action (§. Here Q{t) plays a role of a 
time-dependent potential for the local fermion. In that 
respect the problem in question formally resembles the X- 
ray edge problem solved by Nozieres and De Dominicis 
(ND), with our local fermion being the equivalent to the 
core-hole electron level of the ND problem.l^ This sim¬ 
ilarity has been recognized in the seminal paper by Yu 
and Anderson (YA),I23 in which the authors considered 
the scaling behaviour of a fermionic continuum locally 
coupled to an Einstein phonon, and in which they used 
Hamann’s version of the ND solution.!^ Although our 
model is not fully equivalent to those treated in any of the 
works mentioned above, we can easily adapt their mathe¬ 
matical apparatus to our needs. Despite the fact that the 
resulting effective model is equivalent to the spin-boson 
Hamiltonian,!^ its derivation is much simpler using the 
approach inspired by YA. 

The idea is based on integrating out the local fermion 
field, which is done in imaginary time rather than in the 
energy domain as in To that end we rewrite the full 
system action as (see also Refs. [11],[23]) 

g-S(Q) ^ g-So(Q)^g-A/''drQ(T)<it(r)d(r)^ ^ , 


where the average is taken over the electronic degrees 
of freedom assuming the phonon path Q{t) being fixed. 
The computation of is accomplished in the usual way 
by differentiating it with respect to the coupling constant 
A, perfor ming t he average and integrating again with re¬ 
spect to As a result one obtains 


Z^ = exp 



dTQ{T)D{T, T -I- O"*") 


(9) 


FIG. 5. Phonon displacement probability distribution for 
A = 2.5 and ujq/V = 1 for different inverse temperatnres /3. 


where D{t,t') = — {Trd{T)c[^{ t')) \ is the Matsubara GF 
for the local fermion, calculated in presence of the po¬ 
tential Q{t). It is best computed as a solution of the 
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following exact Dyson equation: 


D{t, t') = Do{t -t') + X [ dr" Dq{t - r") 

Jo 


X Q{t")D{t\t') . 


( 10 ) 


From now on we start making approximations. As was 
realized by ND, in order to access the low-energy be¬ 
haviour of our system it is sufficient to get hold of the 
long-time asymptotics of the solution of Eq. ( [l0| ^ This 
can be conveniently done with the help of the regularized 
version of 


Do{t) 


r 


TTT 



( 11 ) 


where a = [(2/7r) arctan^o]^; = Agmax/E and tq is the 
duration of the hop, which is found from precisely the 
same prescription as the one used by YA. In the strong 
coupling limit the partition function for the system for¬ 
mally equals the one found in m for the Kondo prob¬ 
lem. It possesses a non-trivial scaling behaviour towards 
low energies, which is reflected in a surge of the num¬ 
ber of hopping events. Below the Kondo temperature 
Tx the system undergoes a transition to a new kind of 
ground state, which is universally characterized by Tx- 
From the equivalence of the partition functions we find 
the following relation for the temperature below which 
the bistability vanishes: 

T, « , (16) 


where V denotes the principal value. With this simplifi¬ 
cation the solution of the Dyson equation can be taken 
from [3^. As a result, using the notation ^(r) = AQ(r)/F 
we obtain 


Z^=e^p{- dr[yi(0 + ri(C)] 


( 12 ) 


where 


Di(e) = -(2r/7r)[earctanC- (l/2)ln(l + e2)] (13) 


is an addition to the original harmonic oscillation po¬ 
tential, giving rise to precisely the same terms, which 
are responsible to the adiabatic (Born-Oppenheimer) 
approximation.!^ The other, transient^ term describes re¬ 
tardation effects and represents the next order expansion 
around the adiabatic approximation,!^ 


V 

TiiO = ^ 

TT^ 


dr' 


T — T' 


i + e{ry 


In 


e{r)-e(.r') Ll + e^(^0j’ 


(14) 


The partition function of our system is now a functional 
integral over ^(r) and the oscillator momentum variable 
P{t). It was recognized by Hamann, that in the regime, 
in which the oscillator potential mojQQ'^{T)/2 + Vi{Q{t)) 
develops two distinct minima (bistability regime) the 
dominant paths are given by hopping events between 
them.l^ It was later shown by YA that this picture is 
only insignificantly altered by the kinetic term of the 
oscillator.!^ Following YA and modelling the hopping 
events in the same way one can explicitly write down the 
partition function for the system in question. It turns out 
to be identical to the result in Eq. (61) of YA in which 
the hopping fugacity is replaced by 


y = 


2-a/2 


exp 


—2m 



To 


(15) 


X exp 


-To 




1 

Co 


arctan^o) 


where C is some yet unknown numerical prefactor. When 
adopting the hybridization F as our energy unit Tc be¬ 
comes a single parametric function of x = A^/(FmWg): 
/IcF « (1/Ca:)e’^^/®. From the comparison with the data 
presented in Fig. we conclude that /ScF > 30 and there¬ 
fore we can estimate the constant C as being smaller than 
0.06. 

Bistability signatures have recently been observed 
in nonequilibrium electron transport through molecular 
quantum dots. Although the corresponding simulations 
have been carried out at zero temperature no effects re¬ 
lated to Kondo physics have been seen. It can be at¬ 
tributed to the finite bias voltage playing the role of ef¬ 
fective temperature.!^ On the other hand, there is also 
another reason for absence of Kondo signatures in case 
of small voltages. Starting with an initially decoupled 
system, which is the case in the mentioned nonequlib- 
rium simulations, it takes a finite time tx ~ f/Tx for 
the Kondo effect to become fully established.!^ An es¬ 
timation shows that tc f/Tc is of the same order 
or even larger than timescales addressed by the simu¬ 
lations. Thus the genuine stationary state has not yet 
been achieved, which would explain why the bistability 
signatures were visible. In m the authors reported a 
bistability collapse on a large time scale, which is consis¬ 
tent with our estimation. 

We would like to mention, that although both the nu¬ 
merics as well as the analytical discussion are performed 
for the resonant case ea = 0, we expect everything to hold 
also for finite not too large values of ea. This assertion 
definitely can be shown to hold for the bistability effect 
by a direct computation.!^ As far as the Kondo crossover 
is concerned, it is known to survive in a finite magnetic 
field, which is smaller than Tx- From the mapping be¬ 
tween the models one finds that ea plays the role of the 
magnetic field. Therefore we conclude that for ea < Tc 
one should expect the Kondo crossover to be seen as well. 

Systems of primary interest, in which the bistabil¬ 
ity and YA-type Kondo effect should be observable and 
preferrably also find practical applications are contacted 
molecules. Here the fundamental obstactle is the weak 
electron-phonon coupling. However, this restriction is 
less severe in systems, which are based on carbon nan- 
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FIG. 6 . (Color online) Synthetic quantum system on the basis 
of an ultracold boson-fermion mixture, which is described by 
Q. Two laser fields arranged along the axes ab and cd create 
a deep regular lattice for the fermions, which plays the role 
of the continuum (the electrode in the molecular quantum 
dot realization). Tunnelling between the lattice sites is equal 
and given by 70 . By an auxiliary laser field an impurity site 
is produced, the tunnelling amplitude 7 to which is different 
from 70 . The impurity site is immersed in a BEG cloud, 
with which it can interact and thereby excite lowest lying 
eigenmodes of the cloud. 


otubes. As reported in [ST] the dimensionless electron- 
phonon coupling strength in such systems can exceed 
^ 3, which is sufficient for the bistability as well as Kondo 
crossover to become observable. The only reason why 
that has not yet been observed experimentally is the in¬ 
sufficiently low temperature. This difficulty can certainly 
be overcome in future experiements. 

Another class of systems, in which strong electron- 
phonon coupling is realizable, are ultracold gas mixtures. 
The authors of Ref. [32] have succeded in producing 
a synthetic quantum system with an effective electron- 
phonon coupling strength, which would be sufficient to 
observe the phenomena we have discussed above. In or¬ 
der to see them we envision an optical lattice which can 
be loaded with fermionic atoms, e. g. with ®Li, which 
plays the role of a fermionic continuum, see Fig. Ad¬ 
ditionally a trap for a Bose-Einstein condensate (BEG), 


e. g. of ^^Na atoms, is positioned in the vicinity of one 
of the lattice sites, such that there is at least a partial 
spatial overlap between the BEG in its groundstate and 
one (or several) of the fermionic sites (a prototype sys¬ 
tem is realized in |33|). The latter then plays the role 
of the dot energy level. In presence of boson-fermion in¬ 
teractions there would be a coupling A 7 ^ 0 between the 
lowest-lying harmonic mode of the BEG and the local¬ 
ized fermion level, such that the model Q applies. Due 
to the high tunability of such systems we expect that not 
only the bistability regime might be reached, but also 
the much deeper lying Kondo fixed point. The simplest 
observable of interest is then the spatial dimension of the 
BEG droplet, which is proportional to (Q) and its cumu- 
lants we have discussed above. 

To conclude, we have analyzed the low energy limit of 
the molecular quantum dot model in equilibrium. With 
the help of a path integral Monte Garlo method we have 
computed the probability distribution functions of differ¬ 
ent observables of the harmonic degree of freedom. We 
have shown, that in the strong coupling regime at not too 
low temperature the distribution of the expectation value 
of the average oscillator coordinate becomes bimodal. 
This is a consequence of the bistability effect reported 
previously, which can be described in the framework of 
the Born-Oppenheimer (adiabatic) approximation. Upon 
lowering the temperature the system crosses over into the 
true low energy fixed point, in which the bistability van¬ 
ishes. By means of an analytical expansion around the 
adiabatic approximation we have shown that the emer¬ 
gent behaviour is connected to the Kondo effect, which 
is characterized by a single parameter - the crossover 
temperature - and have related it to the parameters 
of the system under investigation. Furthermore we have 
discussed the experimental implications of the predicted 
phenomena and suggested a setup for their observation 
on the basis of ultracold boson-fermion mixtures. 
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